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CHAPTER I. INTRODUCTION 

For well over a quarter of a ceatury experimentalists (1-8) and 
theoreticians (9-25) have been studying the problem of shock wave 
diffraction, that is, the deflection of a shock wave whose normal path 
has been impeded by some obstacle. Current interest in this problem has 
been generated by researchers (26) investigating the nuclear blast fields 
around aerospace vehicles and around flush-mounted structures (Figure 1) 
in an attempt to accumulate a data.base for survivability and vulnerability 
studies. Such parametric information can be used to determine the nonuni- 
form dynamic loading to be applied in structural analysis programs for the 
design of present day or future generic aerospace systems. 

The interaction of a spherical blast wave with a planar surface, such 
as the examples shown in Figure 1, results in the complete range of shock 
reflections; that is, from regular reflection at 0“ incidence of the 
blast wave with the surface (Figure 2a) to Mach reflection at 90® incidence 
(Figure 2b) . The determination and the understanding of this interaction 
is of importance not only to the structural designer interested in the 
transient blast loading effects but al^ the aerodynamicist interested 
in the mechanics of the flow field . 

The simplest laboratory experiment designed to study the shock 

diffraction problem consists of a two-dimensional wedge or ramp mounted 

on the wall of a shock tube (see Figure 3). Depending on the angle of 

inclination of the ramp with respect to the shock tube wall 0^ and the 

strength of the planar shock (with Mach number M ) , either regular reflec- 

s 

tion or one of the several types of Mach reflection occurs as shoxm in 
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Figure 3. Blast-wave before and after interacting with the compression ramp 
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Figure 3. Regardless of the type of reflection process, this shock 
diffraction problem is self-similar with respect to time since there is 
no characteristic length associated with the problem. 

When the ramp angle is sufficiently large (50“ < 6^ < 90”) regular 
reflection results. As the ramp angle is gradually decreased, the 
shock incident angle (9O°-0^) increases and the regular reflection first 
transitions to a double Mach stem configuration with two triple points 
(see Figure 3) . The second triple point disappears as the ramp angle is 
decreased further and the curvature of the reflected shock reverses. This 
curvature reversal disappears with further decrease in the ramp angle 
and a single Mach stem with a smooth reflected shock appears. For very 
small ramp angles the reflected shock is attached to the ramp edge as 
shown in Figure 3. 

The reason for the formation of a double Mach stem configuration 

during the transition stage from regular to single Mach reflection can be 

explained by a careful examination of the flow field shown in Figures 4a 
* 

and 4b. Let 0^ be the limiting angle for regular reflection. That is^ 
when 0 =0 regular reflection results and when 0 = a tiny 

Mach stem is formed which strikes the ramp perpendicularly. The pressure 
behind the tiny Mach stem (point 3 in Figure 4b) is considerably lower 
than the pressure at point A in the limiting regular reflection case. For 
examples at an incident shock Mach number of 4.71 the pressure at point A 
for the limiting regular reflection case is 147 ^ while the pressure 
behind the tiny Mach stem (point B in Figure 4b) is 62.5. Thus, one Mach 
stem is not sufficient to produce a pressure jump which matches the 


f 



(a) LIMITING REGULAR REFLECTION 



(b) FORMATION OF A DOUBLE MACH STEM 


Figure 4, Transition from regular reflection to double Mach stem 
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limiting reflection value. A second Mach stem is formed such that the 
pressure increase across the second Mach stem matches with the limiting 
regular reflection value. The formation of the second Mach stem is not 
well understood yet* 

The goal of this study is to obtain numerical results for the cases 
of regular reflection and single Mach reflection with a smooth reflected 
shock. The double Mach stem case is not included in the present study. 

All the discontinuities that appear in the flow field are fitted using 
special logic. The reflected shock and the Mach stem are fitted using 
the "sharp shock" technique (27,28). A floating discontinuity fitting 
scheme in conjunction with the method of characteristics is employed to 
fit the slip surface. In the regular reflection as well as the single 
Mach stem case there exists two self-similar stagnation points, that is, 
points at which the self -similar velocity components u-xt and v-y/t are 
zero; the first is located at the juncture of the wall and the ramp 
(saddle singularity), and the second, termed a vortical singularity, is 
located at some point along the ramp (nodal singularity) . In the Mach 
reflection case the slip surface terminates at the vortical singularity 
on the ramp. All streamlines in the self-similar plane converge at the 
nodal singularity, and therefore, the entropy is multivalued. At the 
saddle singularity the streamlines turn away, and the entropy is single- 
valued. The level of entropy on the stagnation streamline and along 
the ramp up to the vortical singularity is equal to that behind the normal 
part of the reflected shock. In the regular reflection case the level of 
entropy between the vortical singularity and the incident shock impingement 


point is equal to that behind the straight part of the reflected shock. 

In the case of single Mach reflection the level of entropy behind the vorti- 
cal singularity and the Mach foot is equal to that behind the Mach foot. 

From an analysis of the equations governing the flow behavior in the 
vicinity of conical, self-similar stagnation points (11,29), it can be 
shown that the pressure is a local maximum at the saddle point of stream- 
lines (juncture of the wall and the ramp) , and this point corresponds 
to a center point of isobars. Similarly it can be shown that the pressure 
is a local minimum at the nodal point of streamlines (vortical singularity) 
which corresponds to a saddle point of isobars. 

In the present study, the two-dimensional, time-dependent Euler 
equations which govern these flows are solved -with initial conditions that 
result in either regular reflection or single Mach reflection of the 
incident shock. The hyperbolic partial differential equations are first 
transformed to include the self-similarity of the problem. Secondly, a 
normalization procedure is incorporated to align the discontinuities as 
computational boundaries to implement the "sharp shock" technique. The 
self-similar transformation reduces these equations from an unsteady to an 
equivalent steady set of mixed elliptic— hyperbolic equations. The 
equations are made totally hyperbolic by reintroducing a time-like or 
residual term which should approach zero in the converged solution. The 
final set of equations is written in strong conservation-law form (30,31) 
and then solved using MacCorraack’s (32) second-order, finite-difference 
algorithm. 
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Unlike previous solutions (12,13,15,22-24) the reflected shock, the 
Mach stem and the slip surface are all treated as sharp discontinuities 
thus resulting in a more accurate description of the inviscld flow field. 
The resulting numerical solutions are compared with available experimental 
data (5) and existing first-order, shock-capturing numerical solutions 
(15,22) . 


CHAPTER II. REGULAR REFLECTION 


When a spherical blast wave strikes a planar surface, regular 
reflection occurs first and then transitions to Mach reflection as the 
shock incident angle increases (Figure 2) . In this chapter regular 
reflection of a planar shock is studied as a prelude to understanding 
more about regular reflection of a spherical incident shock. The simplest 
laboratory experiment designed to study the shock diffraction problem 
consists of a two-dimensional ramp mounted on the wall of a shock-tube 
(Figure 3) . The resulting flow field is self-similar because there is no 
characteristic length associated with the problem. It consists of only 
the reflected shock (Figure 5) which is straight up to the sonic circle 
and then curves to become perpendicular to the shock-tube wall. Between 
the sonic circle and the shock impingement point I the flow field is 
uniform. The flow field linearly grows with time in the physical plane. 

In this problem, there exists two self -similar stagnation points, 
that is, points at which the self-similar velocity components u-x/t and 
v-y/t are zero; the first is located at the juncture of the wall and the 
ramp (saddle singularity), and the second, termed a vortical singularity, 
is located at some point along the ramp (nodal singularity) . All self- 
similar streamlines converge at the nodal point or the vortical singularity, 
and therefore, the entropy is multivalued. At the saddle point the 
streamlines turn away, and the entropy is regular. The ].evel of entropy 
on the stagnation streamline and along the ramp up to the vortical 
singularity is equal to that behind the normal part of the reflected shock, 
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while the level of entropy between the vortical singularity and the 
shock impingement point is equal to that behind the straight part of the 
reflected shock. 

T^ra popular techniques for solving supersonic flow problems are 
currently being used. One is the ’*shock~capturing" method (28), and the 
other is the '’discontinuity-fitting” (27,33) method. The first method does 
not require any special logic to treat the discontinuity and hence yields 
inferior solutions. The discontinuity-fitting procedure requires special 
treatment for all the discontinuities in the flow field (shocks, slip 
surfaces, vortical singularities, etc.). This makes the scheme more 
complicated and involved, but yields a much better solution compared to 
"shock-capturing" results . 

In the present work the "discontinuity-fitting" procedure is adopted 
and the resulting numerical solutions are compared with available 
experimental data (5) and existing first-order, shock-capturing numerical 
solutions (15,22). 


The Transformed Governing Equations 
A Cartesian coordinate system is used in the problem formulation, the 
origin of which is located at the juncture of the wall and the ramp. The 
X— axis is aligned with the wall and the y-axis is normal to the wall and 
in the direction of the ramp (Figure 6a) . Under the assumptions of an 
inviscid, nonheat-conducting, ideal gas, the fluid dynamic equations in 
strong conservation-law form (30,31) for the Independent variable 
transformation t = t, u =n(x,y,t), and ? = C(x,y,t) are 
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where U^E, and F are the conservative variables in the Cartesian 
coordinate system* J is the Jacobian of the transformation. Expressions 
for UjEjF, and J along with the derivation of Equation (1) are presented 
in Appendix A* 

In a shock“fitting procedure the shock is treated as one of the 
computational boundaries, so that jump conditions across the shock can be 
easily applied* This is done through a normalising transformation. For 
the regular reflection problem, the following functions are used for t, 
n and ? which include the self-similarity of the problem and a normali- 
zation of the distance between the ramp and the reflected shock: 


T = t 




X - 3c^(y) 

Xg(y,t) - Xj^(y) 




t 


J 


( 2 ) 


where Xj^(y) represents the equation of the ramp, and x^(y,t) represents 

the equation of the reflected shock. The geometric derivatives 

» 5 , and ? corresponding to the transformation above are used 
t X y t 

in Equation (1). They are derived in Appendix A (see Equations (A13)). 

The self -similar transformation to n and t reduces the unsteady 
gasdynamic equations in the Cartesian system (Equation (Al)), which are 
hyperbolic, to an equivalent steady set of mixed elliptic-hyperbolic 
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equations, that is, in the region between the wall and the sonic circle 
(see Figure 6a) they are elliptic, vdiile above the sonic circle they are 
hyperbolic. The equations are made totally hyperbqlic by reintroducing 
a time-like or residual term (U/J)^ which should approach zero as the 
solution converges. That is, the transformed Euler equations (Equation (1)), 
which are hyperbolic with respect to t, are solved using a time 
asymptotic approach. Because of the self-similarity of the problem, the 
term (U/J)^ approaches zero as x gets large thus establishing a 
convergence criterion. 


Initial Conditions 

The transformation given by Equation (2) results in the computational 

plane shown in Figure 6b. It is bounded by the reflected shock and outer 

boundary, both of which are permeable surfaces, and by the wall and the 

ramp. The region between the wall and the outer boundary is divided into 

^\iax ~ equal Intervals and the region between the ramp and the 

reflected shock is divided into (j - 1) equal intervals . The Inter- 

'’max 

sections of constant-q and constant-C lines generate the discrete 
computational grid used in the finite-difference formulation of the 
problem. Initial conditions (either in terms of flow variables or conser- 
vative variables) are to be specified at all grid points in order to 
initiate the integration of the transformed Equation (1) using MacCormack's 
(32) scheme (Appendix B) . 

To initialize the flow field at time x = 1 given the Incident shock 
Mach number and ramp angle 0^, the pressure and density in region 

(see Figure 6a) are first set equal to unity. The flow conditions in 
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region which are used as the upstream conditions for the reflected 
shock are calculated from the following equations for a moving shock: 


- Cy “ 1) 


^2 Y + 1 


p (Y - 1)M^^ + 2*|pYPiY^^ 

^^2 = 


v, = 0 


P 2 ^ ^^^2 

g = h 


'2 Y “ 3. 


(3) 

(4) 

(5) 

( 6 ) 

(7) 


The subscript 2 in the above equations refers to region (^ • Tlie position 
and the slope of the reflected shock along with the uniform flow conditions 
in region (above the sonic circle in Figure 6a) are then determined 

from the equivalent steady, regular shock reflection equations. This 
procedure is outlined in Appendix C. The conditions in region (^^3^ 
determine the position of the sonic circle at time t = 1. At all points 
along the sonic circle the self-similar velocity is sonic. That is, 

(“3 - + (’'3 - = 

where u^, v^, and a^ are known in the uniform flow region . 

Solving Equation (8) gives the ordinate (x , y ) of any point lying on 
the sonic circle. Nox-?, the outer boundary is chosen such that it falls 
above the sonic circle but below the shock impingement point I (see 


YP_ 

a 2 ^ 

3 P, 


( 8 ) 
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Figure 6a) . Since the outer boundary lies in the uniform flow region 

conditions at grid points along this boundary are exactly known 

using the solution developed in Appendix C. 

The intersection of the sonic circle with the reflected shock is 

determined by simultaneously solving Equation (8) and the slope equation 

for the reflected shock in region . Between this intersection point 

(IP in Figure 6a) and the wall a cubic is used to approximate the 

reflected shock shape. Knowing the shock shape and assuming a self-similar 

flow, that is, Xg = x /t (x is the shock speed), the flow variables 
T ® “t 

behind the reflected shock are given by the following equations, which 
include the Eankine-Hugoniot relations: 


X - x^ 3 

y 
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X = q /cos B 
®t ^ 


( 12 ) 

(13) 

(14) 


B = tan X 


(15) 
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where is the inward normal, is the velocity of the shock in its 

normal direction (see Figure 7) , is the velocity of the flow in the 
direction normal to the shock in region f 2 J , is the shock speed in 

the computational plane, x is the shock speed in the physical plane, 


and Xg is the shock slope. The quantity 


in Equation (16) is 


Vq ~ q. 
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Figure 7. Shock propagation in unsteady flow 
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determined numerically using a second-order central difference formula. 

The subscript s refers to flow conditions behind the reflected shock. 

To initialize the flow field between the ramp and the reflected shockj 
the conditions at the stagnation point (point 0 in Figure 6a) are first 
computed based on the flow conditions behind the normal part of the 
reflected shock (point A in Figure 6a) . 


u = 0 

o 

V = 0 

0 



“a . 


A A 

^V^A 


(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 


Along the ramp between the stagnation point and the sonic circle, a 
parabolic approximation of the flow variables is assumed. The field points 
(l<k<k ) are then initialized by a linear Interpolation 

of the flow variables at the ramp and the reflected shock. Based on the 
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initial flow field, conservative variables U,E,F, and the geometric 

derivatives n. > t] » n » » and £ that are needed in Equation (1), 

t y t X y 

are formulated and stored at each of the grid points. 

Starting with this initialized flow field. Equation (1) is integrated 
(subject to certain boundary conditions discussed in the next section) 
using the explicit, second-order, predictor-corrector MacCormack’s (32) 
scheme. Using a one-dimensional, amplification matrix, stability analysis 
(34) of MacCormack's scheme, a governing integration step size is 
obtained. The integration procedure and step size calculations are 
presented in Appendix B. 


Boundary Conditions 

The computational region is bounded by the reflected shock and the 
outer boundary, both of which are permeable surfaces, and by the wall and 
the ramp, both of which are impermeable surfaces. The boundary condition 
procedures applied at each of these surfaces are discussed below. 

Reflected shock 

The position and the shape of the reflected shock wave are determined 
at each step of the time-asymptotic, integration procedure. The variables 
X , x_ , and x„ which appear in the conservative variables of Equa- 

s Sy St 

tion (1) along with the flow variables at the shock can be determined by 
employing an unsteady variation of the Thomas* "pressure approach*' 

(27,28) for propagating shock waves. In this approach, it is only 
necessary to know the pressure behind the shock in order to alter its 
position for the next time level. The required pressure is obtained by 
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using the normal field point predictor-corrector algorithm at the shock 
but with one-sided differences away from the shock or in the n-direction. 
The shock speed (see Figure 7) and remaining flow variables are given by 
the following equations, which include the Rankine-Hugoniot relations; 
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u “ u„ + 
S 2 


% -^^2 


/ 


1 + 
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(43) 
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1 + x| 


(44) 
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s Y - 1 s 2 


(45) 


The quantity x_ in Equation (34) is determined numerically using a 

second-order central difference formula. The subscript 2 refers to 

flow conditions in the uniform region the subscript s refers to 

flow conditions behind the reflected shock (along k=k in Figure 6a) , 

max ® 

The actual propagation of the shock wave in the numerical procedure 
is accomplished by using a second-order Euler predictor/modified Euler 
corrector 


where 
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The integration step size At is obtained from a stability analysis 
described in Appendix B. 


( 48 ) 
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The above equations are used in the following manner: 

Initially at time step n all flow variables at the shock are known 
including the shock speed and the shape. The pressure behind the shock 
is predicted using the first step of MacCormack^s scheme. The shock wave 
is then moved using Equation (46) . The predicted position permits the 
shock derivatives x to be computed from Equation (34) . The shock 
speedj and other flow variables are then calculated from Equations (32) “(45)* 
The same procedure is followed in the corrector step except that the second 
step of MacCormack^s scheme is used to get the pressure behind the shock 
and Equation (47) is used to correct the shock position. 

Impermeable boundaries 

The impermeable boundaries in the shock diffraction problem consist of 
the wall surface and the ramp surface. Each of these surfaces is aligned 
with a constant coordinate line as a result of the self -similar, normalizing 
transformation. Because of this alignment, and the fact that the flow 
must be tangent to these boundaries, the only variable required at the body 
to advance the field points using Equation (1) is the pressure* However, 
determination of the remaining flow variables and the position of the vorti- 
cal singularity on the ramp is essential in computing the correct surface 
pressure* Discussed below are two different boundary condition procedures 
that were tested for satisfying rhe tangency condition and determining the 
flow variables along the wall and the ramp. 

In the first, a simple Euler predictor/modified Euler corrector with 
one-sided ^-derivatives at the wail and n-derivatives at the ramp for 
Equation (1) is used. The tangency condition itself, that is, v = 0 at 
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the wall and v = u tan 0^ along the ramp, is Imposed after the corrector 
step. Having determined the velocity components from this procedure, the 
self-similar velocities u-x/t and v-y/t are used to locate the vortical 
singularity by noting at what point along the ramp they are identically 
zero. Knowing this location, the appropriate entropy levels are assigned 
to the surface grid points. As mentioned in the Introduction, the level 
of entropy at grid points along the wall and on the ramp up to the 
vortical singularity is equal to that behind the normal part of the 
reflected shock, while the level of entropy at grid points between the 
vortical singularity and the incident shock impingement point is equal that 
behind the straight part of the reflected shock. The corresponding body 
density is obtained from the following expression by using the pressure 
computed by the one-sided finite-difference scheme: 



where s is an appropriate measure of entropy level. The total energy 
e is then recomputed from 


e = 



u^ -f v^ 
2 


(50) 


The second boundary condition procedure tested was that of Kentzer 
(35) . It is based on a method of characteristics approach in combination 
with one-sided finite-differences. Here^ the goal is to derive an 
expression for p^ and valid at the impermeable boundary grid points 
which can be integrated to obtain the surface pressure and the u-component 
of velocity. This procedure is well outlined in Appendix D. Having the 
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u-velocity component, the v-velocity component is computed from the 
surface tangency condition. The self-similar velocities, position of the 
vortical singularity, and the body density are computed in the same way as 
described for the previous boundary condition procedure. 

Using the self-similar property of the flow field in conjunction with 
the surface tangency condition, it can be shown from the normal momentum 
equation that 3p/3n (where n is the direction normal to each surface) 
is zero at the wall and ramp surfaces. Neither Rentzer's scheme nor 
the Euler predictor/modified Euler corrector method satisfy this condition 
exactly because of the approximate one-sided, finite-differences involved. 
Therefore, after the converged solution is obtained using either of the 
above boundary condition procedures, the pressure at the body is recomputed 
after the corrector step to satisfy 3p/9n = 0. This is done in the fol- 
lowing manner. First, the surface normal is drawn and its intersection 
with the first grid line above is found (point R or W in Figure 8) . The 
pressure at the Intersection point (p^ or p^^) is then obtained from a 
simple linear interpolation of the data at two neighboring grid points. 

A simple first-order extrapolation of the form 


satisfies 
and p^ 
stagnation 
pressvires, 


^l,k 


**^,1 % 


I 

1 


(51) 


Sp/3n = 0 to the zeroth order. At the stagnation point both 
are zero. Making use of this condition, the pressure at the 
point is obtained by taking an average of the two extrapolated 
one along the wall and the other along the ramp. 



28 


P 


1.1 



(52) 


A comparison of the different boundary condition schemes is presented 
in the Result section. 


Outer boundary 

The outer boundary (see Figure 6a) is positioned beyond the sonic 
circle (defined by Equation (8)) so that the flow conditions are supersonic 
along it. This allows flow conditions along the outer boundary to be 
specified initially and held fixed during the entire integration procedure. 


Results 

The computational grid for a typical regular reflection case 
consisted of 11 points in the n-direction and 27 points in the ^-direction. 
An average of 300 iteration was required to obtain a converged solution 
and these consumed approximately 15 minutes of computer time on an 
IBM 360/67. 

Numerical results in the form of pressure and density contour plots 
are qualitatively compared with the first-order shock-capturing results 
of Rusanov (15) and Schneyer (22) in Figures 9 and 10, respectively. 
Rusanovas solution was obtained using Godunovas method for an incident 
shock Mach number of 1.89 impinging on a 65° ramp. Most of the contours 
which appear in Figure 9a lie within the captured shock wave, and very few 
describe the flow field between the ramp and reflected shock in comparison 
with the contours of Figure 9b. 


PRESSURE 





(a) RUSANOV (Ref. 15) (b) PRESENT RESULTS 


Figure 9. Comparison of numerically generated pressure contour plots. 

Mg = 1.89, 0J. = 65“ 
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(a) SCHNEYER (Ref. 22) (b) PRESENT RESULTS 


Comparison of numerically generated pressure contour plots. 
Mg - 2, = 63.41® 


3?igure 10. 
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la studying the Mach reflection phenomenon, Schneyer (22) used a 
two-dimensional, Eulerian, hydrodynamic code to obtain the regular 
reflection result shown in Figure 10a. The incident shock Mach number was 
2.0 and the ramp angle was 63.41“. His result exhibits the same quali- 
tative behavior as does Rusanov’s. The present result for the same case 
is shown in Figure 10b. The results of both Schneyer and Rusanov fail 
to reveal the presence of the vortical singularity. 

Law (5) performed a series of experiments on the shock diffraction 
problem for various gases using a Mach-Zehnder interferogram. He tested 
a Mach 4.71 incident shock striking a 60° ramp in oxygen; the result was 
regular reflection. This case in addition to others at the same incident 
shock Mach number but for different ramp angles was obtained numerically 
to demonstrate the flow field behavior in the regular reflection regime. 
The results are presented in Figures 11-16. 

The density and pressure distributions along the wall and the ramp 
are shown in Figure 11. At the stagnation point (point C of Figure 11), 
the density and pressure reach a local minimum, while at the vortical 
singularity (point D of Figure 11) , the pressure is continuous and at a 
local minimum, and the density is discontinuous . A partial plot (see 
Figure 12) of the self-similar velocity along the ramp reveals the two 
self-similar stagnation points at A (juncture of the wall and the ramp) 
and B (vortical singularity) . 

In Figure 13, results from the different body boundary condition 
procedures are compared. Both the Euler predictor /modified Euler 
corrector and Kentzer’s scheme yield very nearly the same results. The 




Figure 12, Self-similar surface velocity distribution. Mg = 4.71, 
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Figure L3. Comparison of ramp surface pressure distribution for different boundary condition 

procedures, llg = 4.71, 0^ = 60° 




Figure 15. Self-similar velocity directional plot. Mg = 4.71, 0 
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Figure 16. Comparison of experimental and computed shock positions. M = 4.71, 0 
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oscillations near the stagnation point are a result of the one-sided, 
finite-differences used in these schemes. Imposing 3p/3n - 0 seems to 
yield a much better solution without any oscillations in the flow variables 
near the stagnation point. 

Pressure and density contour plots of these computational region are 
shown in Figure 14* The centerpoint of isobars near the wall-ramp inter- 
section point, and the saddle point of isobars near the vortical 
singularity (for which moving away from the vortical singularity the 
pressure increases along the ramp and decreases perpendicular to the ramp) 
can be clearly observed in the figure. In the density contour plot, the 
convergence of the various isopycnics at the vortical singularity can be 
observed. The behavior of the flow near the stagnation points in this 
unsteady two-dimensional self-similar problem exhibits the same behavior 
as the steady, self-similar, three-dimensional flow about an external 
axial corner (36) • 

The self-similar streamline pattern can be visualized by observing 
the velocity vector directional plot of the computational plane shown in 
Figure 15. Notice that all the streamline converge at the vortical 
singularity . 

A comparison of the interferogram obtained by Law (5) with the 
numerically computed shock shape is shown in Figure 16. If an overlay of 
the two results were made by matching shock impingement points, the 
experimental shock location would fall inside the numerical solution. 

The reason for the discrepancy is probably twofold: First, the viscous 

effects (the majority of which can be observed near the wall-ramp 
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intersection point) might have the effect of decreasing the ramp angle as 
a result of the boundary layer growth with distance from the shock 
impingement point. The reduced ramp angle in turn results in a smaller 
shock standoff distance. Second, the computed solution assumes flow of an 
ideal gas (y = 1.4). Thus, high temperature effects on the internal 
energy such as molecular, vibrational excitation are not taken into account. 

The effect of varying the ramp angle for a given shock Mach number 
of 4.71 on the shock standoff distance (r ), position of the vortical 

SO 

singularity (r ) , location of the sonic circle (r ) , and shock impinge- 
ment point (r^) are shown in Figure 17. The standoff distance exhibits 
almost a linear variation with ramp angle between the limit for regular 
reflection and the last computed case of 0^ = 85°. The vortical singtxlar- 
ity moves towards the wall with increasing ramp angle and actually attaches 
itself to the wall for values of 0^ greater than 77°. The location of 
the sonic circle along the ramp, and the shock impingement point are 
identical at the limit for regular reflection. As 0^ increases 
the sonic circle moves toward the wall while the impingement point moves 


away from the wall. 



Figure 17, Variations of geometrical variables with ramp angle. M. 
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CHAPTER III. SINGLE MACH EEELECTIOH 

A typical single I&ch reflection of an incident shock is shown in 

Figure 18. The self-similar flow field is somewhat complicated in this 

case by the existence of a triple point at which the reflected shock, the 

Mach stem and the incident shock meet. Emanating from the triple point 

is a slip surface which intersects the ramp at the vortical singularity. 

In Figure 18, M denotes the self-similar Mach number. A sonic line 
ss 

exists in most of the single Mach reflection cases in the region between 
the reflected shock and the slip surface. Below this sonic line (region I 
in Figure 18) and in the region between the Mach stem and the slip surface 
(region III) the self-similar Mach number is subsonic (M < 1), while 
above the sonic line (region II) it is supersonic (M > 1) . 

In ^h^^g problem, there are two self -similar stagnation points, that 
is, points at which the self-similar velocity components u-x/t and v-y/t 
are zero; the first is located at the juncture of the wall and the ramp 
(saddle point), and the second, termed a vortical singularity, is located 
at the point where the slip surface meets the ramp (nodal singularity) . 

All the self-similar streamlines converge at the vortical singularity and 
the entropy is multivalued. v^aiuc of entropy on the stagnation 

streamline and along the ramp up to the vortical singularity is equal to 
that behind the normal part of the reflected shock, while the entropy 
between the vortical singularity and the Mach foot is the same as that 
behind the foot of the Mach stem. 

In the present work, the reflected shock and the Mach stem (including 
the triple point) are fitted using the "sharp shock" (27,33) technique. 



WALL 


Figure 18. Single Mach reflection of the incident planar shock 
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A floating discontinuity-fitting scheme in conjunction with the method of 
characteristics is developed to fit the slip surface. 


Double Normalization Procedure 

A Cartesian coordinate system is used in the problem formulation with 
the origin located at the juncture of the X'/all and the ramp. The x-axis is 

aligned with the wall and the y-axis is normal to the wall (Figure 19a) . 

The gasdynamic equations in this Cartesian system are given by Equation (Al) 
for the assumptions stated in Appendix A. 

In order to apply the "sharp shock" technique, the reflected shock and 
the Mach stem are used as computational boundaries. This is done by means 
of a double normalizing transformation. The following functions are used 
for T, ri, and 5 which include the self-similarity of the problem, a normal- 
ization of the distance between the ramp and the reflected shock and a 
normalization of the distance between the wall and the Mach stem: 


T = t 




X - 

" Kg(?,T) - 

= 1 


(53) 


where Xj^(S,t) represents the equation of the ramp, represents 

the equation of the reflected shock and Y^(n,T) represents the equation of 
the Mach stem. Note that the ramp and the shock shapes are defined in 
terms of the computational variables n, snd x and not in terms of 
X, y, and t. Prom Figure 19a it can be seen that such a representation 
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is necessary because constant n and constant | lines are not parallel 

to the Cartesian x and y axis, respectively. 

The equations corresponding to the independent variables t, n, and 

I are given by Equation (1) in strong conservation-law form. The geometric 

derivatives n , n > n » j £^nd ? appearing in this equation are 

t 3C y ti y 

derived in Appendix A (see Equations (A15) and (A16)). 

The self -similar transformation to n and ? reduces the unsteady 
gasdynamic equations in x> y> and t s^^stem* which are hyperbolic, to an 
equivalent steady set of mixed elliptic-hyperbolic equations; they are 
elliptic in regions of subsonic self-similar velocity 0^^^ 1) ^tid 

hyperbolic in regions of supersonic self -similar velocity (M > 1) . The 

5S 

equations are made totally hyperbolic by reintroducing a time-like or 

residual term ((U/J) ). Because of the self-similar nature of the flow 

X 

field this tlme-like term should approach zero in the converged solution. 

Initial Conditions 

The transformation given by Equation (53) results in the computational 

plane shown in Figure 19b. It is bounded by the reflected shock and the 

Mach stem, both of which are permeable boundaries, and by the wall and the 

ramp. The coordinate ris Is zero at the ramp and equal to one at the 

reflected shock. Similarly K Is zero at the wall and equal to one at 

the Mach stem. The slip surface floats within the mesh generated by the 

double normalization (Figure 19a) . The region between the ramp and the 

reflected shock is divided into (j - 1) equal intervals and the 

max 

region between the wall and the Mach stem is divided into 

equal intervals. This determines the mesh spacings An and A^. The 
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intersections of constant n and constant C lines generate the discrete 
computational grid used in the finite-difference formulation of the problem. 
Initial conditions are specified at all grid points in order to initiate 
the. integration of the transformed Equation (1) using MacCormaclc^s (32) 
scheme (see Appendix B) • 

To initialize the flow field at time t = 1 given the incident shock 
Mach number M^ and the ramp angle 6^, the pressure and the density in 
region © (see Figure 20) are first set equal to unity. The flow condi- 
tions in region © , which are used as the upstream conditions for the 
reflected shock are calculated from Equations (3) to (7) . Referring to 
Figure 20, an initial value for the triple point trajectory angle (x) 
is assumed. Corresponding to this assumed value of x> the triple point 
solution is computed from an equivalent steady approach described in 
Appendix E. This gives the flow conditions at points 3 and 4 lying on 
either side of the slip surface at the triple point (see Figure 20) • The 
reflected shock slope the Mach stem slope and the slip surface 

angle a are also obtained from the triple point solution. 

Assuming some standoff distance for the reflected shock (distance 
0-A) at time t = 1, a cubic is used to approximate the reflected shock 
shape between the triple point and the wall. This cubic satisfies the 
conditions that the shock he normal to the x^all at point A and the slope 
at the triple point be equal to that determined by the triple point 


solution (tan ij)„) . Similarly, assuming some value for the distance between 
the origin 0 and the Mach foot B, a cubic is used to approximate the Mach 
stem between the triple point and the ramp. This cubic satisfies the 



Figure 20^ Floi^ field initialization 
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conditions th-'Jt the Mach stem be normal to the ramp at the foot and the 
slope at the triple point be equal to that determined by the triple point 
solution (tan . The slip surface is initially approximated by a 
straight line with a slope tan (90®-a) . This straight slip surface meets 
the ramp at the point denoted by VS in Figure 20. 

Even though the double nomnalization requires that both the reflected 
shock and the Mach stem be represented in terms of the computational 
variables n, 5, and t the calculation of the flow variables behind them 
and their actual propagation requires a representation in terms of the 
physical variables x, y, and t. The reflected shock is represented by 

X = X (C,t) = X [y(n=l,^),t] (54) 

where x (y,t) is the representation in terms of the physical variables. 
Similarly 3 the ifech stem is represented by 

y = x^(n,T) = yg[x(c=l,Ti),t] 

where y (x, t) is the representation in terms of the physical variables. 

s 

Knowing the initial reflected shock shape and assuming a self-similar 
floWj that is, 

X^(C,t) 

Xs (e,T) = (55) 

T 

the flow variables behind the reflected shock are given by the following 
equations, which include the Rankine-Hugonio t relations (see Figure 21) : 
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Figure 21. Reflected shock propagation in unsteady flow 
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where n is the inward reflected shock normal, q is the velocity of 
s s 

the shock in its normal direction, is the velocity of the flow normal 

to the shock in region , X is the shock speed in the computational 

plane, x is the shock speed in the physical plane, and x is the 

s 

t y 

shock slope. The quantity X (C,t) appearing in Equation (61) is 

computed numerically using a second-order central difference formula. 

The subscript 3s refers to flow conditions behind the reflected shock 

(along n = 1 in Figure 19b) . 

The flow conditions behind the Mach stem are computed in a similar 
fashion knowing the initial shape and assuming a self-similar flow 
(see Figure 22). 
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Figure 22. Mach stem propagation in unsteady flow 
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where n is the outward Mach stem normal, q is the velocity of the 
s s 

Mach stem in its normal direction, Y is the shock speed in the 

s 

T 

computational plane, y is the shock speed in the physical plane, and 

y is the Mach stem slope. The quantity appearing in Equation (76) is 
s 

X 

computed numerically using a central difference formula. The subscript 1 
refers to flow conditions in region , and the subscript 4s refers to 
flow conditions behind the Mach stem (along C 1 in. Figure 19b) . 

With the flow conditions^ behind the reflected shock and the Mach stem 
known, all the field points (l<k<k ) ate now initiali 3 ed. 

The conditions at the stagnation point (point 0 in Figure 20) are computed 
from Equations (25) to (31) based on the normal part of the reflected 
shock (point A) . The point where the slip surface meets the ramp is a 
vortical singularity where the self-similar velocity components are zero. 
Thus, initially the velocities at point VS (see Figure 20) are assumed to be 
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(87) 

V = — ^ — I 

vs T J 

Along the ramp between point 0 and point VS and between point VS and 
point B (Mach foot) the velocities are linearly interpolated and the 
pressure is approximated by a parabola. Along the wall (A— 0) a parabolic 
approximation of the flow variables is used. The entropy along the wall 
and the ramp up to the p^-int VS is equal to that behind the normal part 
of the reflected shock, and between VS and the Mach foot it is equal to 
that behind the foot of the Mach stem. With pressure and entropy kno\m 
along the ramp and the wall, the density is computed from 



P 



( 88 ) 


where s is an appropriate measure of the entropy level. The total 

energy e is then computed. The pressure at the field points (1 < k < k^^^, 

1 < 1 < i ) are obtained by a lijrear interpolation of the pressure at 

"’max 

the reflected shock and the ramp. The pressure along the slip surface is 
then obtained by a linear interpolation using the values at the neighboring 
grid points. The side of the slip surface facing the reflected shock is 
denoted by "a" (see Figure 20) and the side facing the Mach stem is 
denoted by "b." The pressure on either side of the slip surface is the 
same but the velocities are not. The velocity components along the slip 
surface on side "a" are obtained by linear interpolation using the values 
at point VS and point 3 at the triple point. Similarly the velocity 


components along the slip surface on side "b" are obtained by linear 
interpolation using the values at point VS and point 4 at the triple 
point. The u-velocity component along the slip surface on side "a" is 
recomputed to satisfy the jump condition 

(u^ - Uj^) tanOO” - a) = (89) 

The entropy along the slip surface on side "a" and side "b" is equal to 
Sg and respectively (the slip surface is a self-similar streamline along 
which the entropy is constant) . Knowing the entropy and the pressure, the 
density along the slip surface is computed. Based on the flow conditions 
along the slip surface on side "a” and the reflected shock, the field 
points lying in region I (Figure 20) are initialized. Similarly, based 
on the flow conditions along the slip surface on side "b" and the Mach 
stem, the field points lying in region II are initialized using linear 
inter p o lation . 

Starting with this initialized flow field. Equation (1) is integrated 
(subject to certain boundary conditions discussed in the next section) 
using the explicit, second-order, MacCormack's (32) scheme. Since the slip 
surface floats within the n,5 mesh system a floating- fitting scheme in 
conjunction with the method of characteristics is developed to propagate 
the slip surface. Under this scheme differencing across the slip surface 
is forbidden. Thus, special one-sided differencing formulas (37) are used 
at grid points neighboring the slip surface. The floating-fitting scheme 
along with the special differencing formulas are explained in a later 


section 
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Boundary Conditions 

The computational region is bounded by the reflected shock, the Mach 
stem, the wall and the ramp. The boundary condition procedures used at 
each of these surfaces are discussed below. 


Reflected shock 

The position, shape and the speed of the reflected shock wave are 
determined at each step of the time-asjrmptotic integration procedure. 

The variables X X (Cjt), and X (C,t) which appear in Equation (1) 

S 5 — S 

5 T 

along with the flow variables behind the shock are determined by employing 
the unsteady version of the Thomas' "pressure approach" (27,28) for 
propagating shock waves. As mentioned in Chapter II, in this approach 
it is necessary to know only the pressure behind the reflected shock 
(p^g) in order to alter its position for the next time level. This required 
pressure is obtained by using the normal field point predictor-corrector 
algorithm at the reflected shock but with one-sided differences in the 
n-direction . As mentioned in the previous section (Initial Conditions) 
in order to compute the shock speed X (C,t), it is necessary to define an 

t 

equivalent reflected shock shape in terms of the physical variables 
t, X, and y. Such a representation is given by Equation (54). Knowing 
pressure the remaining flow variables are given by the following equations, 
which include the Rankine-Hugoniot relations. 
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The velocity components u and v and the total energy e are then 

oS <dS oS 

computed from Equations (68) to (70) • The actual propagation of the 
reflected shock along with the Mach stem is presented under a separate 
subheading . 


Mach stem 

The variables Y (n,T)j Y (n^x), and Y (n,x) along with the flow 
s s s 

n T 

variables behind the Mach stem are determined from pressure (p ) using 

^■S 

the same Thomas' "pressure approach" employed for the reflected shock. 

The pressure behind the Mach stem (p ) is obtained from the finite- 

HS 

difference algorithm using one-sided differences in the g-direction. The 
remaining flow variables are given by the following equations: 
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The velocity components and and the total energy are 

computed from Equations (84) to (86) . 

Impermeable boundaries 

The boundary condition procedure used at the wall and the ramp is 
exactly the same as that used for the regular reflection case except that 
Kentzer's scheme (35) was not used because it gives the same results as 
one-sided finite-differences. 

Shock Speed Calculations 

The actual propagation of the reflected shock and the Mach stem 
in the numerical procedure is accomplished by using a second-order Euler 
predictor/modified Euler corrector. For the reflected shock it is given by 

x“'*'^(5,t:) = x/(C,t) (?,t)At (104) 

s s s 

T 

Xg’*'^?,T) - Xg^(C,T) +1 (C,T) + x"‘^kC,T)^AT 


(105) 
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For the Mach stem 
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Equations (104) and (106) are the predictor step and Equations (105) and 

(107) are the corrector step. It is necessary to represent the shocks in 

terms of the physical variables in order to evaluate the shock speeds 

X (C,t) and Y (n,T). They are evaluated in the following manner; 
s s 

T T 


X3 (I.T) - Z (y.t) + i y 
Tty 


(108) 


n=i 


Yg (n>T) = y^ (x,t) + y^ 

T t X 


c=i 


(109) 


where 



(5,t) 

J6 





(110) 


( 111 ) 


= 5Y^ (i,t) 
^ n=i T 


( 112 ) 


= (1,t) + (l-n)X^ (1,t) 


(113) 


'c=l 
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(1,t) = tan(90“ - (0,t) 


X 


= Xg(l,T) -Xj^( 1 ,t) 


(114) 

(115) 


C“1 


= Y^(1,t) 


(116) 


n=i 


X and y are evaluated from Equations (60) and (75) , respectively. 

^t ^t 

It can be seen from Equation (108) that evaluation of X (?,x) requires 

y which in turn requires Y (1,t) . Similarly, evaluation of 

■^fn=i 

Y (n»x) requires X (1 ,t) which is the reflected shock speed at the 
s s 

T T 

triple point. Since the triple point moves with a Mach number M , X (l,x) 

3 S 

T 

is simply the speed of the incident shock wave* 


X^ (l,x) = M /T 


(117) 


Substituting n ~ 1 in Equation (109) an expression for Y (1,t) is 

I 

obtained 


^ ^s ^s 

T t X X 


(118) 


In Equation (118) y and y are evaluated at the triple point. 

s . s 

t X 

This completes the calculation of shock speeds. Knowing X_ (^,x) and 

®x 

Y (HjT) the reflected and the Mach stem are advanced using Equations (104) 

S 

T 

to (107). 
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Floating-Fitting Procedure for the Slip Surface 

In a usual discontinuity-fitting procedure, the discontinuity is 
transformed into a computational boundary by means of a normalizing 
transformation. According to Moretti (37) this is not necessary. His 
floating-fitting procedure allows one to float the discontinuity within 
the existing mesh and still fit it by using special one-sided differencing 
at grid points neighboring the discontinuity. The idea is not to allow 
differencing across the discontinuity. The actual propagation of the 
discontinuity in Moretti’ s approach is done using Kentzer’s scheme. In 
the present analysis, Moretti ’s floating— fitting idea is used to treat 
the slip surface. Instead of using Kentzer’s scheme, the method of 
characteristics is used to compute the flow conditions along the slip 
surface on either side (38) . 

Figure 19a shows the slip surface floating within the existing mesh 
system. Since differencing across the discontinuity either in space or 
time is strictly forbidden in a fitting approach, special one-side differ- 
encing formulas are used at grid points neighboring the discontinuity, 
instead of the usual equally spaced difference approximations (MacCormack’s 
scheme uses forward differences in the predictor and backward differences 
in the corrector). Application of the usual MacCormack’s scheme at the 
grid point 0 neighboring the slip surface in Figure 23 requires conservative 
variables at grid points 5 and 6 lying on the other side of the slip 
surface. Since this is forbidden, the forward, differencing in the 

n-direction (E in Equation (AlO)) is modified to the following (Refer- 

n 

ence 37) : 
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(119) 


SE 


t^^SE- - (4 - 5.^)E^ + (1 - e^)E^1 


2At) 


( 120 ) 


SE 


, ' 


1 - e 


S = 


1 1 + e 


( 121 ) 


2 n 1 


(122) 


Equation (119) requires two backward grid points (1 and 2) . If only one 
backward grid point is available then Equation (120) is modified to 


E 


SE 


" ^ [(^)® ^ r 


fr- ®SE - 


] 


(123) 


If no backward grid point is available then Equation (119) is modified to 

(124) 


^SE “ S 


e An 

n 

Similarly, the backward differencing in the ^-direction (F^ in Equa- 
tion (AlO)) at grid point 0 is modified to 

[Fsx - (1 - - ^2^3^ 




AC 


Vc 


sx 


where 


2AfI 


(125) 


( 126 ) 
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6 


1 


1 - e. 


1 + e. 


(127) 


6^ = (128) 

Equation (125) requires two forward grid points (3 and 4) . If only one 
fortrard grid point is available then Equation (126) is modified to 


f 


C 


SX 


AC 



1 + e, 




(129) 


If no forward grid point is available then Equation (125) is modified to 


F 


C 


0 



(130) 


Similar special differencing formulas are used at grid points neighboring 
the slip surface on the other side. 

All of these special one-sided formulas require the evaluation of 
the conservative variables at points where the slip surface intersects the 
constant n and constant C lines. In order to formulate the conservative 
variables, the flow conditions along the slip surface must be evaluated at 
each time level. This is done using the method of characteristics. 

The flow field is initially assumed on either side of the slip surface 
at points where it intersects the constant n and constant C lines. As 
integration proceeds in the time (x) direction, the location of the initially 
assumed slip surface keeps changing along with the flow variables on 
either side until the correct self-similar solution is reached. The 
actual propagation of the slip surface is carried out only along the 
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constant n lines, using the method of characteristics. The location and 
the flow conditions at points where the slip surface intersects the con- 
stant C lines are then obtained by a linear interpolation of the values at 
two neighboring slip surface points lying on constant n lines (see 
Figure 24). In Figure 24, points and "bj" represent two sides of 

the slip surface at a constant n line at the initial time level "n." At 
the new time level (t + At) n+1, they are given by "a" and "b.” The problem 
here is to locate this new slip surface position and to compute the flow 

conditions at "a" and "b." Out of the ten flow variables (p , p , u , v , 

a a. a a 

®a* ^b’ ‘^b’ “b’ ^b ®b^ ^^a’ ’^a* ^b’ “b* ^b^ 

evaluated. The densities p and p, and the total energies e and e, 

a D a o 

can be obtained from 


P 


a 



pb 



(131) 


(132) 


e 

a 



P_(u ^ 
a a 


+ V 


Pb , + 'b^> 

7 ^-" 2 


(133) 


(134) 


where s and s, are some measure of the entropy values on either side 
a b 

of the slip surface. They are the same as the values at the triple point 
because the slip surface is a self-similar streamline along which the 


entropy is constant. 
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STREAMLINE 
CHARACTERISTICS 
DOWN RUNNING 
CHARACTERISTICS 
fTl 

CONSTANT 97 -LINE 


COMPUTATIONAL 
GRID POINTS 


CONSTANT ^-LINES 



CONSTANT 
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Figure 24. Method of characteristics for the unsteady slip surface propagation 
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Evaluation of the pressures (p and p, ) and the velocities 

& P 

(u , u, , V , and v, ) require six simultaneous algebraic equations. As 
aba’ b 

previously noted, the slip surface is advanced only along the constant 
n lines. Thus, the characteristic compatibility relations are derived 
only in the (C - t) plane. Figure 24 shows the slip surface location at 
the old (n) and new (n+1) time level at a constant n line. From the new 
time level location, the C and C characteristics are drawn which strike 
the old time level at points 1 and 2, respectively. The compatibility 
relation along the C characteristics is given by 


(Pa ■ 


pcC 




(u - u.) + 


peg. 


r 


+ e. 


(V^ - v^) 


up + pc^nu + pc^n V + 

n ^ n y n 


pc?, 


X 


pc? 

+ 2- 




(Y? + uv ) 
V p n/ 


/e ^ + V ^ 

/ X y 


At 


/n P \ 

1 + uu 1 

V P n/ 


(135) 


'X y 

Similarly, the compatibility relation along the C characteristics is 
given by 




pc? 


<Pt - Pj) - 




i '"b - “z’ ■ 




2 ^ ('b - -z’ 
2 + ^2 


up + pc2q u + pc2p V “ _ 

X n y /i 


r 2 + £ 2 \ p n/ 

’x y 
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/ S 
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(136) 
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The line joining the old slip surface location (point a^ in Figure 24) 
and the new slip surface location (point a) is nothing but the streamline 
characteristics in the (5 - t) plane. Compatibility relations are derived 
along the streamline characteristics on either side of the slip surface. 


On side "a” it takes the form 


£(u ~u )— S(v - V ) = 
a ar a al^ 




(137) 


The similar relation on side "b" is 


■ ’'b 


i' ■ ■ ['y('^ * ““n) ■ 


—'•h uv ) At 


(138) 


The jump conditions across a moving slip surface are 


Pa = Pb 


(139) 


q •n.=q. *n,=q 

si b si ^sl 


(140) 


where q^ and q^ are the velocity vectors on either side of the slip 

surface, n , is the slip surface normal and q , is the velocity of the 
si ^sl 

slip surface in its normal direction. The slip surface is defined by 


y - = 0 


(141) 


Using Equation (141) , the second jump condition is rewritten in the form 


where fyg-j^) is the slope of the slip surface 


(142) 


\ 
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The unhnotTO flow variables p^, p^, u^, v^, and are now 

obtained by solving Equations (135) -(139) and (142) simultaneously. The 
actual propagation of the slip surface is carried out by 


where 




4 -^ 


is the speed of the slip surface and is evaluated in a 


manner similar to the shock speed. When the flow field converges, the 

slip surface speed (y ^ ) should converge to y /t. 

\ ®1/t ®1 


Results 

The computational grid for a typical single Mach reflection case 
consisted of 6 points in the n-direction and 31 points in the i; -direction. 
All average of 400 iterations was required to obtain a converged solution 
and required approximately 15 minutes of computer time on an IBM 360/67. 

Humerical results in the form of pressure contours are qualitatively 
compared with the first-order shock-capturing results of Rusanov's 
solution in Figure 25. Rusanov's solution was obtained using Godunov's 
method for an incident shock Mach number of 1.89 impinging on a 30“ ramp. 
Moat of the contours which appear in Figure 25a lie within the captured 
shock waves, and very few describe the flow field bounded by the reflected 
shock, the Mach stem, the wall and the ramp in comparison with the contours 
of Figure 25b. 

Law (5) performed a series of experiments on the shock diffraction 
problem for various gases using a Mach-Zehnder interferometer. He tested 
two cases which resulted in single Mach reflection. The ramp angle for 
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both the cases was 40° and the incident shock Mach number for one case was 
1,89 and for the other 2,10, Numerical results were generated for these 
two cases to demonstrate the flow field behavior in the single Mach reflec- 
tion regime. The numerical results are presented in Figures 26-35, 

The density and pressure distributions along the wall and the ramp for 
two cases are shown in Figures 26 and 27. The juncture of the wall and 
the ramp is a stagnation point (point C) at which pressure and density 
reach a local maximum. The point where the slip surface meets the ramp is 
a vortical singularity (point D) at which the pressure is continuous and 
reaches a local minimum. The vortical singularity is nothing but a slip 
surface at a point at which the density takes a jump because of the 
discontinuous behavior of the entropy. The numerical results clearly 
exhibit this flow field behavior as predicted by Ludlof f and Friedman (11) . 

The mesh in the physical plane is automatically generated fay the 
double normalizing transformation. As reflected shock and the Mach stem 
change their shapes during the iteration process, the mesh in the physical 
plane also keeps deforming until the self-similar flow field is established. 
Figures (28) and (29) show the converged mesh in the pfavsical plane for 
incident shock Mach numbers 1.89 and 2.1, respectively, for a ramp angle 
of 40°. The slip surface is clearly seen to float within the physical 
mesh. 

Pressure contour plots of the physical region are shown for two 
cases in Figures 30 and 31. The centerpoint of Isobars near the wall-ramp 
juncture, and the saddle point of isobars near the vortical singularity can 
be clearly observed. By doing a local analysis of the gasdynamic equations 





Figure 28. 


Mesh in physical plane. Mg = 1.89, 6j- = 40^ 



Figure 29^ Mesh in physical plane. Mg = 2,1^ ~ 



Figure 30* Pressure contours* Mg =1*89, ©r ” 
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Figure 31* Pressure contours. M - 2*1, 0- = 40 






Figure 34. Comparison of experimental and computed shock and slip surface positions. 
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Figure 35. Comparison of experimental and computed shock and slip surface positions. 
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at these singularities (stagnation point is a saddle singularity and the 
vortical singularity is a nodal singularity) Ludloff and Friedman (11) 
came up with the same behavior for isobars as seen in the numerical 
solutions. Figures 30 and 31 also exhibit the continuous behavior of 
pressure across the slip surface. In order to show any discontinuous 
behavior of the flow field as a sharp jump in the contour plot, the 
contour program requires that such a discontinuity be treated as one of 
the boundaries of the computational region because of the various inter- 
polations involved. Since the slip surface is floated within the 
computational mesh the contour program cannot bring out the true sharp 
jump in the density across the slip surface in a density contour plot. 

The density contour plot might look as though the slip surface was captured 
within a mesh interval. 

The self-similar velocity directional plot for two cases are shown 
in Figures 32 and 33. The self-similar streamline pattern can be visualized 
from these plots. Notice that all the streamlines tend to converge at 
the vortical singularity (nodal singularity). The streamlines also diverge 
away from the stagnation point (saddle point) . Only the stagnation 
streamline passes through the stagnation point. 

The comparison of the interferogram obtained by Law (5) with the 
numerically computed shock and slip surface shape is shown for two cases 
in Figures 34 and 35. The triple point trajectory angle (x) in the numeri- 
cal solution is larger than that shown in the experimental interferogram. 

The reason for the discrepancy is probably two fold: First, the viscous 

effects (the majority of which can be observed near the wall-ramp 
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intersection) might have the effect of decreasing the ramp angle as a 
result of the boundary layer growth with distance from the Mach foot. The 
reduced ramp angle in turn results in larger triple point trajectory angle. 
Second, the computed solution assumes flow of an ideal gas (y = 1.4). 

Thus, high temperature effects on the internal energy such as molecular, 
vibrational excitations are not taken into account. The slip surface in 
the numerical solution comes out to he nearly straight as seen in the 
experimental picture. In addition, a small self-similar supersonic region 
lies between the slip surface and the reflected shock. The sonic line 
bounding this supersonic region is shown in the numerical results in 
Figures 34 and 35. 




CHAPTER IV. CONCLUDING REMARKS 


The discontinuity-fitting procedure developed in this report for com- 
puting the shock diffraction problem for the regular and the single Mach 
reflection is capable of accurately predicting the inviscid flow field with 
its reflected shock, the Mach stem, the slip surface and the vortical singu- 
larity. The solution in the neighborhood of the self-similar stagnation 
points exhibit gasdynamic equations with regards to the behavior of the 
self-similar streamlines, isobars and isopycnics. The present numerical 
results are a considerable improvement over the early first-order numerical 
solutions and compare favorably with available experimental data. 

The present work treats only the regular reflection and the single 
Mach reflection cases. In order to develop a discontinuity-fitting proce- 
dure for the double Mach reflection case, a good a priori knowledge of the 
flow structure is required. Thus, a development of a good shock-capturing 
solution for the double Mach stem case is very desirable to understand 
what exactly is going on. 

Extension of the present planar shock diffraction problem t<' the 
spherical shock diffraction problem is suggested. 
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APPENDIX A. STRONG CONSERVATION-LAW FORM OF THE GOVERNING 
EQUATIONS AND THE GEOMETRIC DERIVATIVES 


Since the equations of motion in fluid mechanics are derived from con- 
servation principles (mass, momentum, and energy), it is often convenient 
to cast the equations in divergence form or conservation-law form which 
explicitly displays the conserved quantities such as mass, momentum, and 
energy. In the Cartesian system (x,y,t) the gas-dynamic equations (contin- 
uity, x-monentum, y-momentum, and energy) for invlscid, uonheat-conducting, 
and adiabatic flow can be written in conservation- law form as : 


where 


Dt + E^ + Fy = 0 


(Al) 



u and V are the velocity components in x and y directions, and p, p, and 
e are the pressure, density, and total energy per unit volume. The system 
of equations is made complete by specifying the total energy in the form; 

e = y ^ ^ P ^ ^ 

In most fluid mechanics problems it is necessary to make a coordinate 
transformation from the Cartesian coordinates (x, y, t) to some other system 
(n, t), in order to facilitate the easy application of surface boundary 

conditions on arbitrary shaped bodies* Sometimes coordinate transformations 
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are needed to incorporate some of the special features of the flow field 
(conical flows, self-similar flows, etc.) in the numerical formulation of 
the problem. 

Let the new coordinates t, t1i and ^ be related to the Cartesian sys- 
tem t, X, and y by the transformation: 

T = t 

t| = n(x, y, t) y (A3) 

? = ?(x, y, t) 

To arrive at the transformed equations, the derivatives with respect to x, 
y, and t in Equation (Al) are replaced in terms of the derivatives with 
respect to t, n, and C in the following manner: 
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^ _a_ , X jr 

3x 3n 'it 




-i- = -L + X t 
3y 3n ’^y 3c ^y 


X = X + X 

3t 3 t ^ 3n 




3C 


Making use of Equation (A4) , Equation (Al) can be written as 

ft + itUn + + n**n * Vc + V" * V? * ° 

Equation (A5) can always be rewritten as 


+ E^ + + H' =0 


where 


U’ = U 

E’ = HfcU + + Hyf 

^t^ 
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5” = t E + EF 
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+ Eti„ + Fn„ + UCt + + Ff 
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(AS) 


(A6) 


(A7) 
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All the terms in the untransformed Equation (Al) are derivatives of the 

unknown f our-component vectors (U, E, and F) with respect to the independent 

variables (t,x,y). This is said to be in strong conservation-law form (30). 

The transformed Equation (A6) is said to be in weak conservation-law form 

because of the presence of an undifferentiated term . This term is 

analogous to the fictitious body force term* The presence of the H' term 

in the governing transformed equation is undesirable for two reasons* First, 

it prevents the achievement of overall conservation of mass, momentum, and 

energy. Second, it involves several second derivatives ^ r\^ , r\^ ^ y 

^T) 5 

5x s £v ) * analytical expressions required to evaluate these second 

derivatives may be difficult to obtain. As a result these are evaluated 
numerically thus increasing the computer time. 

In the present work, the transformed Equation (A5) is rewritten i 
strong conservation-law form to avoid the undesirable features of the weak 
conservation-law form. In order to bring Equation (A6) into a strong 
conservation-law form the H’ term must somehow be removed by including 
appropriate terms into U’ , and F’ before the derivative is taken. 

The technique of Viviand (31) is applied here. 

The Jacobian of the transformation is given by 

hx Htl 


9 (n > € a *r) „ c 
3(x,y,t) 


^x^y ■" ^y^x 


V '^t 


If the transformation is regular the Jacobian is neither zero nor infinite. 
Assuming the Jacobian to be finite, all the terms in Equation (A5) are 
divided by the Jacobian J. It can then be rearranged as 
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L Jt u Jt, l_ 


+ E? + FC 
N 




+ E' 

It can be easily shown that all the terms inside the fourth bracket cancel 
out. Thusj Equation (A9) is composed of only the first three bracketed 
terms. The strong conservation-law of the transformed equations can thus 
be written in a simpler fashion as: 


+ E„ + Tj = 0 


(ilO) 


where 




(All) 


U = U/J 

E = (UtIj. + EtIj^ + Friy)/J 
F = (U5^ + E?^ + F?y)/J 

^ " ’^y^x 

For an analytical transformation the geometric derivatives n^j Hy, E}-> 
and Cy can be evaluated analytically. For a numerical transformation 
these geometric derivatives will have to be evaluated numerically. 

In the regular reflection problem (refer Figure 5) the independent 
variable transformation t = t, p = p(x,y,t), and g = 5(x,y,t) which 
includes the self- similarity of the problem and a normalization of the dis- 
tance between the ramp and the reflected shock is given by: 

T = t 


n = 


Xg(y,t) - x^(y) 


^ t 




(A12) 
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where S|^Cy) = y cot 0^ represents the equation of the ramp, and Xg(y,t) 
represents the equation of the reflected shock. The geometric derivatives 
required by Equation (AlO) are: 




hf- = 

s b 
1 






X. 




” ’^(^s - ) 

V y y/ 






^x = 0 


«y-T 


(A13) 


Since = 0, the Jacobian reduces to J = hj^Sy. 

In the Mach reflection problem (refer Figure 19a) , the transformation 
involves a double normalization procedure. The transformation functions 
T, ri) and C include the self-similarity of the problem, a normalization of 
the distance between the ramp and the reflected shock and a normalization 
of the distance between the wall and the Mach stem. They are given by; 


T = t 


X - 


n = 


C = 


XjjCe.T) 


XgC?,T) - x^(c,t) 

y - Yt,(x) 


(A14) 


YgCn,x) - y^(t) 

Since the wall is aligned with the x-axis^ the equation of the wall is just 
Y^(t) =0 or C 0* In Equation (A14) , == y cot 0^ - ^yg(0,T)cot 0^ 

represents the equation of the ramp (n “ 0) , represents the equation 

of the reflected shock, and YgCn^x) represents the equation of the Mach 
stem. The body and the shock shapes are defined in terms of the computational 
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variables n, C, and t and not in terras of the physical variables x, y, 
and t. Such a representation is necessary because the constant ri and 
constant 5 lines are not parallel to the x and y axis, respectively. 


Corresponding to the transformation given by Equations (A14), the 
geometric derivatives are obtained as follows: 




h,. = 




3(n,Y,t) 

3Cn,g,T) 

9(x,y,t) 


n„ = 


^x 




3(x,n,t) 

9(n»g»T) 

3(x,v,t) 

9CTi,g,'i:) 

3(x,y,n) 

3(n,g,T) 

3(x,v,t) 

3Cn,g,T) 

3(S»y. t) 
3(n,g,x) 
3(x,y,t) 
3(n,g,T) 

3(x,g,t) 

3(n,g,x) 

3(x,y,t) 

3(n,g,T) 

3(x,y,g) 

3(n,g,T) 

3(x,y,t) 

3Cn,g,x) 


yg ^x 


t 

g X 






X Vj^ - x.y 


X V- - X-V 


x.y - X y_ 
g-'x T~^g 

X y^ - x,.y 

n g g n 


Vg ■ Vn 


X 

T1 

X y^ - x^y 


X y - X y 
T^n n T 

X y. - x»y 




(A15) 


t 
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Evaltiation of (A15) requires the following: 

34 = nx ci,t) + (1 - n)x, ( 5 ,x) 

T T 

= hXg (C,t) + Cl - (C,t) 

= XgC?,T) - X^(C,T) 

= SYg (q,x) 

T 

Vjr = Yg(n,T) 

Yn = 

n 


(AI 6 ) 


/ 
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APPEiroiX B. INTEGRATION PROCEDURE AND STEP SIZE CALCULATIONS 


MacCormack (32) has constructed a two step, preferential, predictor- 
corrector sequence for use in solving systems of differential equations 
written in the conservation-law form. The scheme is second order in both 
time and space* In application to nonlinear equations with several depen- 
dent and independent variables, the method has low storage requirements and 
simple programming logic* 

As applied to Equation (AlO) MacCormack^ s method is as follows: 

(Bl) 


0IX+1 -n\ 

^j,k An\^Vi.k j.k/ A?\^j,k+l 

"j.k - 2 [’’j.k + “j.k - -ST (®J.k - "j-l.kj - 4?i^J.k - 'j.k-ljJ 


(B2) 


The tilde that appears over certain of the variables denotes the predicted 
value of that particular variable* The subscripts j and k refer to mesh 
indices whereas the subscript n refers to the time* 

In this version forward differences are used in the predictor and back- 
ward differences in the corrector. However, one could use backward differ- 
ences in the predictor and forward differences in the corrector* Another 
possibility is to use a forward difference for the n-derivative and a back- 
ward difference for the ^-derivative in the predictor and the opposite in 
the corrector. Because of these various options MacCormack^ s scheme is 
termed a preferential difference scheme . 

In the case of a boundary mesh point for which the forward grid is not 


available, the forward difference in the predictor in that direction is 
modified to a backward difference in that direction. Similarly, if the 
backward grid is not available then the backward difference in the corrector 
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is modified to a forward difference in that direction. For example, for 

all the grid points along the ramp no backward grid point is available in 

the n-direction. At these grid points the term E, , - E, . in the 

corrector is modified to . 

3^1, k 3,k 

The integration step size At must be specified to initiate the calcu- 
lation. The maximum allowable step size Ax^ in the n-direction and the 
maximum alloxmble step size Ax^ in the ^-direction are obtained from the 
one-dimensional, amplification matrix, stability analysis (34) of MacGormack 
scheme. They are given by 



= CN 

An 

(B3) 

1 '^max, n j 

At^ 

= CN 

AC 

(B4) 

U .1 


max , 5 1 


where CN is the Courant number, a is the maximum eigenvalue in the 

max,n 

(n-t) plane, and a r maximum eigenvalue in the (5-x) plane. 

max, ^ 

For the calculation to be stable, the minimum of the two step sizes Ax 


and At^ is used: 


Ax = minCAx^, Ax^) 


(B5) 


In order to compute these maximum eigenvalues first the equations of 
motion are written in nonconservation form in terms of the transformed 
coordinates variables r\, and x. The continuity equation and the energy 
or the entropy equation are coupled together to eliminate any derivatives of 
density. This is done in the following manner: 


Continuity: 

+ q * Vp -f pV • q = 0 

(B6) 

Energy: 

pj. - c^p^ + q ■ (Vp - c^Vp) =? 0 

(B7) 


f 
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where 


q = 

V ^ 


oil + VJ 


9 * 3 

ax ^ 3y 


c2 = lE. 

P 

Multiplying Equation (B6) by and adding it to Equation (B7) results in 

+ up^ + vpy + pc2(u^ + Vy) = 0 (B8) 

In terms of the transformed coordinate variables Pj and t Equation (B8) 
becomes 

p^ + Pj^u + p^v + pc2(u^^^ + + v^Py) = 0 (B9) 

where 

u = n + up + vp (Bio) 

t. X y 

V = + ug^ + vgy (Bll) 

The x-momentum and the y-momentum equations are also written in terms of 
the transformed coordinates. 

x-momentum: u_ + p p /p + p„E /p + u 13 + u_v = 0 (B12) 

T *^p X “^g^x p g 

y-momentum; v_ + p p /p + p_g /p + v u + v^v “ 0 (B13) 

^ n y C y n E 

Equations (B9), (B12) , and (B13) are written in matrix form as: 

Qt + ^iQti ^ 2 % = 0 


where 



(B15) 



101 


The matrix A, has three eigenvalues the maximum of which is a 

1 max, 71 

Similarly the matrix A 2 has three eigenvalues the maximum of which is 
o . The eigenvalues of are obtained by solving the matrix equation 

[Aj - lal = 0 (B16) 

Where I is the identity matrix. Solving Equation (B16) yields the follow- 
ing three eigenvalues: 

Ai 

rr -L 


U 


— / ^ 

*^2,3 = “ - V 'lx + 'Ij 


The absolute maximum is given by 


max,n 


= |u| + + r\y^ 


(B17) 

(B18) 

(B19) 


Similarly solving [a, - Icr| =0 yields the following three eigenvalues: 


A' 

= V 


A2 


^3 = v±c/i7n; 


and 


W.cl ■ hi + V + «y^ 

The integration step size is now given by 

An 


At = minlCN 


u + c/nj^ + n,^ 


CM 




-^1 + + r 


(B20) 

(B21) 

(B22) 

(B23) 


'X 'y ■ ’ V -X "y 

Equation (B23) is evaluated at each of the grid points in the computational 
plane and the smallest value of At over all the grid points is then 
chosen as the integration step size. The Courant number CN is usually 
chosen to be one or slightly less than one. 
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APPENDIX C. EXACT SOLUTION FOR REGULAR REFLECTION 

When a planar shock strikes a wall it will reflect in one of two forms, 
regular reflection or Mach reflection. The form that occurs depends on the 
shock strength and the shock incident angle. In the present problem (refer 
Figure Cl) the incident planar blast wave denoted by its strength Mg 
strikes the ramp with an incident angle of (90-flj.), where 0^ is the ramp 
angle measured from the positive x-axis. For incident shock Mach numbers 
greater than 1.5, regular reflection results as long as the incident angle 
is less than 39“ (39). 

In the numerical formulation of the regular reflection problem, the 
computational region is chosen such that the outer boundary (refer Figure 6a) 
falls between the sonic circle and the point I where the incident shock 
strikes the ramp. Along the outer boundary exact t\?o-dimensional regular 
reflection results are specified and kept fixed throughout the iterative 
process. The exact regular reflection results are obtained by making use of 
various shock relations in the folloxiring manner. 

As the incident shock moves with a Mach number Mg, the shock incident 
point I (Figure Cl) moves up the ramp with a Mach number M^/cos 0j-. The 
shock relations such as the Rankine Hugoniot jump conditions are applicable 
only when the shock is at rest. These shock relations can be applied to a 
moving shock by merely employing a moving coordinate system relative to 
which the shock is at rest. By placing a moving coordinate (x',y‘) rigidly 
attached to the moving point I, the stationary region in the (x,y) 
system becomes nonstationary in the (x',y’) system. With respect to the 
moving system (x*,y') the nonstationary flow in region is parallel to 




I 
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the ramp and has a Mach number M^/cos 0j. as shown in Figure C2. The 
transformation from (x,y) to Cx’jy') system alters only the velocities in 
the regions^^^^,^^^^^, The pressure and density remain unchanged. The 
velocities with respect to the Cx',y') system are denoted by a prime. 

To obtain the flow variables in regions and^^^^, first the pressure 

and density in region^^^^l^ are chosen to be unity (i-e-i P^ “ Pj == !)• Then 
the following equations found in NACA 1135 (40) are used. 


A- 

~ = v'f speed of sound in region/ IJ 

= M_/y ~ velocity of the incident shock 

S 5 


1 cos 0, 




velocity in region ( 1 Jwith 
to aystem ^ 


respect 


?2 = Pi 


P2 " Pi 


'2YM„^ - (y - 


Y + 


(Y + l)Mg^ 

(Y - 1)M^2 4 . 2 


]■ 


pressure in region 


© 

density in region 


(Cl) 

(C2) 

(C3) 

(C4) 

(C5) 

(C6) 


0 , = f - 9^ 


Mj,^ = M[ sin - M, 


, - 1) 3.) _ 

'^7 (y + dXiMI" 


6 = tan 


“1 


2 cot 0 j(m/ ~ 1) 


L2 + MJ2(Y + 1 “ 2 sin2 0^) J 


(C7) 

(C8) 


velocity in region (2 J with (C9) 
respect to (x'jy') system 


flow deflection angle (CIO) 

from region to 




Figure C2. Regular reflection in a moving x’,y' 
Cartesian system rigidly attached at the point I 
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~ speed of sound in region 


ion(|^ 


cell) 


(C12) 


In region the flow again becomes parallel to the ramp. Thus knowing 
6 and the reflected shock angle Sg is found by solving the following 
polynomial; 

sin® 0„ + b sin*^ 0„ + c sin^ 0„ + d = 0 (C13) 


where 


+ 2 

^ _ _ Y 3^2 5 


m; 


1 2 


(C14) 


2M^^ + 1 

c + 

1*2 


d = - 


cos^ d 

MJ" 


il±12i+ V - 3. 


M, 


t2 


sin^ 6 


(C15) 


(C16) 


Equation (C15) has three roots, the smallest of which corresponds to a 
decrease in entropy and should therefore be disregarded according to the 
second law of thermodynamics. The largest root corresponds to the strong 
shock. The middle root, which corresponds to the weak shock, is the one of 
interest. 

^ sin 0g (C17) 


= *12 


?3 *= T>2 


1 - 


4(lf g - 1) (YM|g + 1) 


(Y + l)2M|gM^2 


~ velocity in region 0 
respect to (x’y ) syste 


with 
system 


- (Y - 1) 


Y + 1 


pressure in region 


0 


(CIS) 


(C19) 
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(.y + 


(Y - 1)M§2 + 2 


density in region 


© 


(C2Q) 


Knowing the velocities q*, and q^ in the (x%y^) system, the velocities 
q^j q^, and q^ in the (x,y) system are obtained by employing the simple 
trans formation : 


where 


- 

. qj - S . Ujl + Vj3 

(C21) 

->■ -5-r 

<12 = ^2 - 

• = u^i + v^3 

(022) 

<13 “ 'Is - 

. A A. 

• 'll “ “3^ ^3^ 

(023) 

= q^(-cos 

0J.X “ sin 0j.|) 

(024) 

= q^(~cos[0j. + 6]i - sin[0^ + 6]j) 

(C25) 

= q^(-cos 

0j.i - sin 9^5) 

(026) 


and Cu^,v^)j and (u^jV^) are the Cartesian velocity components in 

regions ©•© , and respectively. 
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APPENDIX D. KEKTZER’S SCHEME FOR IMPERMEABLE BOUNDARIES 

An Impermeable boundary is one across which no mass can flow such as a 
solid surface or a plane of symmetry. At an impermeable boundary a surface 
tangency condition must be satisfied. Proper implementation of the surface 
boundary condition is a crucial step in computing the correct body pressure 
distribution. One method of applying the surface boundary condition is to 
use Kentzer's (35) scheme at the body grid points. 

Kentzer’s scheme is based on the method of characteristics approach in 
combination with one-sided fin±te differences. Here, the aim is to derive 
an expression for p^ valid at the body points which can be integrated in 
a predictor-corrector fashion to obtain the body pressure. This is achieved 
by combining the characteristic compatibility relation and the surface 
tangency condition in differential form. The procedure is outlined below 
for both the ramp and the wall. 

The eigenvalues of the time dependent Euler equations have already been 

derived in Appendix B (see Equations (B6) through Equation (B21))- The left 
Ai 

eigenvectors y^ corresponding to the eigenvalues of the A^ matrix are 
obtained by solving 

y^HA^ - lojl) =0 i = 1, 2, 3 (Dl) 

Similarly solving 

y^2(^2 - Ia^2) =0 i = 1, 2, 3 (D2) 

yields the left eigenvectors y^^ corresponding to the A 2 matrix. The 


final result is 
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At 

= (0, Tiyj - Tljj) 




- 1 1 . ± 


pen 


X 


pen. 


A 


+ n„^ 




yf = (0. Cy, -5^) 


p = 

' 2,3 


X 

pc?, 


= 1 , ± 


X 


pc?^ 


/? ^ + ? 

/ X 5 


/? 2 + ? 2 
/ ^x y 


(D3) 




The compatibility relations are obtained from these eigenvectors and 
eigenvalues . 

Referring to Figure Dl, only the down running characteristics drawn in 
the C? - t) plane strikes the wall grid point. The eigenvalue associated 
with this characteristic is Kentzer's scheme requires only the com- 

patibility relation along this down running characteristic. The compatibil- 
ity relation is derived by starting from Equation (B14) . 


Qt + = 0 


(D4) 


Multiplying Equation (D4) throughout by y^^ and making use of Equation (D2) 
it results in the form 




(D5) 


Substituting for y^^ from Equation (E3) and for Q from Equation (B15) , 


Equation (D5) simplifies to 


I 


^ Z + ^ Z ^ /g2 + g2 ^ 

X y V ^x y 


= - up^ + ^x^n ■*■ ^y^r\ 


'”=^K /V„ . . 

- - ■ . ’ — ( 4- uu. 


2 ^ ^ 2\ P 


c ^ + c 

X ^y 




Along the wall the surface tangency condition in differential form is given 


v^ = 0 


V = 0 

n 


In addition along the wall (plane of S3nmnetry) v is zero and is zero. 

Combining Equation (D6) and Equation (D7) and then substituting for u^ 
from Equation (B12) yields the following expression for p^ valid only at 
the wall grid points. 


/a? a? , _ . 9 

~ + P^Vn p — 


X - 
— - uu„ 
T p n 


Equations (D8) and (D9) are integrated in a predictor^corrector fashion to 
get the pressure and the u-velocity at the wall grid points at the new time 
level. 


~n+l n , , ^n . 

"l.k - "j,k <PT>j.k 


-n+l n , / A 


predictor 


(DIO) 
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n+1 n . 1 r, ..n -- 

^j,k “ Pj,fc 2 

n+i n . 1 1 / j_ /•" 

“■J,k = “j.k+2 ['Vj,k+ K)j.kJ 


> corrector 


(Dll) 


In evaluating and forward differences are used for the n and C 
derivatives in the predictor. In the corrector backv-xard differences are 
used for the n derivatives and forward differences for the g derivatives* 

Knowing the pressure and the u-component of the velocity all the other 
flow variables can be easily computed. This procedure is outlined in 
Chapter II tinder Boundary Conditions. 

Similar to the analysis presented above* an exrpression for p^ and 
are now derived for the ramp grid points. Referring to Figure E2* only the 
doT\m running characteristics drawn in the (n - t) plane strikes the ramp 
grid point* The eigenvalue associated with this characteristics is cr^^. 

The compatibility relation along this down running characteristics is 
obtained by multiplying Equation (D4) by y^^. 


" -yt^A,Q 


3 


(D12) 


Ai 

Substituting for from Equation (D3) and for Q from Equation (B15), 


Equation (D12) results in 

pcn^ 




A 


-t" n 

^ y 




pen. 


-[ 


vpg + " ~J 


JA + "i 

pen 


(v^ + a^lv^) 


X 


+ n. 


2 ^ + ™c) 


_ p^, Y 


(D13) 


Ill 


Along the ramp the surface tangency condition in differential form is given 
by 

v^ = tan 0- 


(D14) 


tan 0^ 


In addition, u is zero along the ramp. Combining Equation (D14) and then 
substituting for u^ from Equation (B12) yields the following expression 
for p^ valid only at the ramp points. 




At 

0_^p 

3 


pen. 


A 




At 

- 

3 n 


pen^ 


A 


+ n 2 




pen 


+ pc^e u_ + pc^e v_ - 


X 




y c 


/n ^ + n ^ 

V X y 


n - 


c 

p 


pen 


2 y p 

/ n ^ -f n ^ 

/ X y 


■) 


^ = _p p 


(D15) 


(D16) 


Equations (D15) and (D16) are integrated in a predictor-corrector fashion 
described by Equations (DIO) and (Dll) . In evaluating and f on^rard 
differences are used for n and g derivatives in the predictor. In the 
corrector backward differences are used for the g derivatives and fonvard 
differences for the n derivatives* Knowing the pressure and the 
u-component of the velocity all the flow variables are easily obtained. 

By combining the compatibility relation with the surface tangency con- 
dition in differential form, the disadvantages of the true method of charac- 
teristics, the iterations and the interpolations to get the data at specific 
points on a characteristic, are eliminated in Kentzer*s scheme. 


DOWN RUNNING 
CHARACTERISTIC 
IN ^“T PLANE 



CONSTANT 
PLANE 


WALL" n. j,k 

Figure Dl. Kentzer's scheme at the wall point (k = 1) 


RAMP 



Figure D2. Kentzer's scheme at the ramp point (j = 1) 
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APPENDIX E. EXACT TRIPLE POINT SOLUTION 


t-Ihea the incident angle (that is the angle between the incident planar 
shock and the ramp) is greater than 39°, Mach reflection occurs as long as 
the incident shock Mach number is greater than 1.5 (39). As mentioned in 
the Introduction (Chapter I) , the Mach reflection can take various forms 
depending on the incident angle and the incident shock strength. The pres- 
ent problem considers only the single Mach reflection case in which only one 
triple point is present. Since the flow field is self-similar the triple 
point moves along a straight line denoted by the triple point trajectory 
angle x Figure El. 

As was pointed out in Appendix C, all the shock jump conditions are 
true only if the shock is at rest. In order to obtain a solution to the 
moving triple point where the incident shock, the reflected shock, the Mach 
stem, and the slip surface meet, a moving Cartesian coordinate (x’,y') is 
rigidly placed at the moving triple point. With respect to the (x',y*) 
system the triple point is at rest and the flow comes into the triple point 
along the triple point trajectory with a Mach number Mg/cos(6j. + x)> as 
shown in Figure E2. The transformation from (x,y) to (x',y*) system alters 
only the velocities in regions and^^^^. The pressure and 
density remain unchanged. The velocities in the two system (x,y and x',y') 
are related to each other by means of a simple transformation. 

The triple point solution is first obtained in the x*,y’ system where 
all the shock relations found in NACA 1135 (AO) are applicable. Knowing flow 
variables in region region © is easily obtained using oblique shock 

relations. To solve regions uniquely , an iterative procedure is 
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necessary to satisfy the two jump conditions across the slip surface. In 


the x^,y^ system the slip surface is at rest. The two jump conditions 
across a stationary slip surface are one, the pressure must be same on either 
side of the slip surface and two, the velocity vector on either side must be 


parallel to the slip surface. Knowing region 


Lon^^l^ , 


the following procedure 


explains how to obtain the flow variables in regions ^^^2^ , and ©• 

Pj " P]^ = 1 ~ pressure and density in region 

= j ~ ~ speed of sound in region - (E2) 

q = M„/Y^ ~ velocity of the incident shock (E3) 

s ® 

ql = y - .^ , r ~ velocity in region f©^with respect (E4) 

1 cosCo-r + XJ-/-i i\ Ivy 
t to Cx ,y ) system ^ 


MJ = — 

X a^ 


r2YM^ - (T - 1)1 
- - 

r (T + i)m| -] 
^^L(y - 1)m| + 2J 


pressure in region 


density in region 


agion(^^^^ 

ion 


0 = - V 

1 2 ^ ^ 




4(m2j - i)(ym|i + 1) 

(Y + l)^MiiM{2 

2 cot 01 (Mg^ - 1) 

2 + M[2(y -1-1-2 sin2 0^) 


(ElO) 


(Ell) 


(E12) 


^2 

m: = — 


(E13) 


This completes region calculations . The following iterative procedure 
determines region region 0 flow conditions* 

1* For a given Macti number there is a maximum flow deflection angle, 
Knowing in region^^^^^ the maximum flow deflection angle 

the reflected shock is computed from 




^2^ ^max - ^ 




tan 6, 


(E14) 


. X T L irt " 


(E15) 


The flow deflection angle 5^^ across the reflected shock has to be less 
than or equal to 5^^^. 

2. An initial value for the Mach stem angle is chosen. The Mach 

stem being a strong shock, the initial value for Cjj is chosen to be 89.99* 

3. Corresponding to the assumed Mach stem angle and the Mach 

number M|, the floxi? deflection angle is obtained from 


2 cot %(MJ2 sin^ - 1) 

6w = tan ^ — 

2 + MJ2(y + 1 - 2 sin2 

4, The pressure in region ^4^ is computed from 

Ptt = Pl|l +V^ (Ml^ % - 1) j 


(E16) 


(E17) 
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5. The slip surface angle a is given by 


a “ 90° - Sjj - 0,. - X 


(E18) 


6 . Since the flow in region has to be parallel to the slip surface 

the flow deflection angle across the reflected shock is given by 

5r = ct - (90° - 0,. - X - (E19) 

7 . If is greater than given by Equation (D14) , the initial 

guess for is reduced by 0 . 01 ° and the calculation is repeated from 

step 3 until 6 ^^ becomes equal to or less than 

8 . Knowing the flow deflection angle and the Mach number M^ the 

shock angle is computed. The procedure is outlined in Appendix C. 


9. The pressure in region^^^S^is then computed from 


P, = P. 


; 1 Sin^ - 1) 


(E20) 


10. Across the slip surface the pressure must be same (i.e., Pg = p^^) . 
If the convergence criteria 

|p 3 “ P 4 I < 10"^^ (E21) 

is not satisfied then the assumed value of Cjj is reduced by 0.005 and the 
calculation is repeated from step 3. This repetition is continued until the 
convergence criteria satisfied. 

11. The total velocities in region and q^ in region are 

given by 


qi = q’ 


< = 'll 


1 - 



4(M^^ sin2 - 1)(YM^2 ^ Tj 


*2 


(Y + 1)^M2**^ sin^ 5 


4(M|^ sin^ - 1) (YM{^ sin^ Cjyi + 1) 
(Y + sin2 


(E22) 


(E23) 
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12. Knowing the velocities q’, q^, q*, and qj^ in x*,y' system, 
the velocities q^, qg, q^, and q^ in x,y system are obtained by employ- 


ing the simple transformation; 

^2 “ “2^ 

^3 = ^3 - = “ 3 ^ + V t^26) 

^4 = - q{ = “ 4 ^ + (E27) 

where 

qj = q|[-cos(0j. + x)i - sin(0^ + x)j] (E28) 

^2 ~ *l2l“cos(0j. + X + 5)i - sin(0^ + X + fi)j] (E29) 

^3 ~ ®1^ (E30) 

94 “ <j4^“®^^ (E31) 


and (UjjVj^), (u 2 ,V 2 ), (ugjVg), and (uj^,v^) are the Cartesian velocity 
components in regions 

13, The slopes of the reflected shock and the Mach stem are given by 


'J’M " ®r + X % 

(E32) 


(E33) 


2\f 3\ andf 4\ respectively. 
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Figure El* Mach reflection in a fixed 
x,y Cartesian system 


Figure E2* Mach reflection in a moving x' ,y^ 
system rigidly attached at the triple point TP 
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